home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / NRPAS13.ARJ / SOLVDE.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  3KB  |  98 lines
  1. PROCEDURE solvde(itmax: integer; conv,slowc: real; scalv: glscalv;
  2.        indexv: glindex; ne,nb,m: integer; VAR y: glyarray;
  3.        nyj,nyk: integer; VAR c: glcarray; nci,ncj,nck: integer;
  4.        VAR s: glsarray; nsi,nsj: integer);
  5. (* Programs using routine SOLVDE must define the types
  6. TYPE
  7.    glindex = ARRAY [1..nyj] OF integer;
  8.    glscalv = ARRAY [1..nyj] OF real;
  9.    glyarray = ARRAY [1..nyj,1..nyk] OF real;
  10.    glcarray = ARRAY [1..nci,1..ncj,1..nck] OF real;
  11.    glsarray = ARRAY [1..nsi,1..nsj] OF real;
  12. in the main routine. *)
  13. LABEL 99;
  14. CONST
  15.    nmax=10;
  16. VAR
  17.    err,errj,fac,vmax,vz: real;
  18.    ic1,ic2,ic3,ic4,it: integer;
  19.    j,j1,j2,j3,j4,j5,j6,j7,j8,j9: integer;
  20.    jc1,jcf,jv,k,k1,k2,km,kp,nvars: integer;
  21.    ermax: ARRAY [1..nmax] OF real;
  22.    kmax: ARRAY [1..nmax] OF integer;
  23. BEGIN
  24.    k1 := 1;
  25.    k2 := m;
  26.    nvars := ne*m;
  27.    j1 := 1;
  28.    j2 := nb;
  29.    j3 := nb+1;
  30.    j4 := ne;
  31.    j5 := j4+j1;
  32.    j6 := j4+j2;
  33.    j7 := j4+j3;
  34.    j8 := j4+j4;
  35.    j9 := j8+j1;
  36.    ic1 := 1;
  37.    ic2 := ne-nb;
  38.    ic3 := ic2+1;
  39.    ic4 := ne;
  40.    jc1 := 1;
  41.    jcf := ic3;
  42.    FOR it := 1 TO itmax DO BEGIN
  43.       k := k1;
  44.       difeq(k,k1,k2,j9,ic3,ic4,indexv,ne,s,nsi,nsj,y,nyj,nyk);
  45.       pinvs(ic3,ic4,j5,j9,jc1,k1,c,nci,ncj,nck,s,nsi,nsj);
  46.       FOR k := k1+1 TO k2 DO BEGIN
  47.          kp := k-1;
  48.          difeq(k,k1,k2,j9,ic1,ic4,indexv,ne,s,nsi,nsj,y,nyj,nyk);
  49.          red(ic1,ic4,j1,j2,j3,j4,j9,ic3,jc1,jcf,kp,c,nci,ncj,nck,s,nsi,nsj);
  50.          pinvs(ic1,ic4,j3,j9,jc1,k,c,nci,ncj,nck,s,nsi,nsj)
  51.       END;
  52.       k := k2+1;
  53.       difeq(k,k1,k2,j9,ic1,ic2,indexv,ne,
  54.          s,nsi,nsj,y,nyj,nyk);
  55.       red(ic1,ic2,j5,j6,j7,j8,j9,ic3,jc1,jcf,k2,
  56.          c,nci,ncj,nck,s,nsi,nsj);
  57.       pinvs(ic1,ic2,j7,j9,jcf,k2+1,
  58.          c,nci,ncj,nck,s,nsi,nsj);
  59.       bksub(ne,nb,jcf,k1,k2,c,nci,ncj,nck);
  60.       err := 0.0;
  61.       FOR j := 1 TO ne DO BEGIN
  62.          jv := indexv[j];
  63.          errj := 0.0;
  64.          km := 0;
  65.          vmax := 0.0;
  66.          FOR k := k1 TO k2 DO BEGIN
  67.             vz := abs(c[j,1,k]);
  68.             IF (vz > vmax)  THEN BEGIN
  69.                 vmax := vz;
  70.                 km := k
  71.             END;
  72.             errj := errj+vz
  73.          END;
  74.          err := err+errj/scalv[jv];
  75.          ermax[j] := c[j,1,km]/scalv[jv];
  76.          kmax[j] := km
  77.       END;
  78.       err := err/nvars;
  79.       fac := 1.0;
  80.       IF (err > slowc) THEN fac := slowc/err;
  81.       FOR jv := 1 TO ne DO BEGIN
  82.          j := indexv[jv];
  83.          FOR k := k1 TO k2 DO BEGIN
  84.             y[j,k] := y[j,k]-fac*c[jv,1,k]
  85.          END
  86.       END;
  87.       writeln;
  88.       writeln('Iter.':8,'Error':9,'FAC':9);
  89.       writeln(it:6,err:12:6,fac:11:6);
  90.       writeln('Var.':8,'Kmax':8,'Max. Error':14);
  91.       FOR j := 1 TO ne DO writeln(indexv[j]:6,
  92.          kmax[j]:9,ermax[j]:14:6);
  93.       IF (err < conv) THEN GOTO 99
  94.    END;
  95.    writeln('pause in routine SOLVDE');
  96.    writeln('too many iterations'); readln;
  97. 99:   END;
  98.